###################################################################################################

#Figure3_PanelB
#This figure plots the evolution of the average value of commodities rceived from January to November of 2017 by treatment and control
###################################################################################################

### Setup #######################################################

#DEFAULT COLORS
fColor <- "black"
mColor <- "gray60"
lightColor <- "#97B8DE"
darkColor <- "#4669BD"


ration.value_total <- c("totalrec","totalnonrec")
ration.name_total  <- c("total reconciled","total non-reconciled")


month.list <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun","Jul", "Aug","Sep","Oct", "Nov")

rect_r <- data.frame(xmin=7, xmax=10, 
                     ymin=-Inf, ymax=Inf)

fontSize <- 16
theme_set(theme_gray(base_size = fontSize))


#######################################################################################
# Reconciled
#######################################################################################
tmp.outcome <- "value"
seqalong <- seq_along(ration.value_total)

#Gathering data for reconciled
r <- 1
tmp.ration <- ration.value_total[r]
tmp.ration1 <- ration.name_total[r]

#treatment

tmp.data.c <- readstata13::read.dta13(paste(SurveyDataDir,tmp.ration,"_",tmp.outcome,"_meanreceipts_reconciliation_time_series_treatment0.dta",sep = ""))
tmp.data.t <- readstata13::read.dta13(paste(SurveyDataDir,tmp.ration,"_",tmp.outcome,"_meanreceipts_reconciliation_time_series_treatment1.dta",sep = ""))


#reshape the dataset to include values for month and whether predicted/observed
meltdf.c <- reshape2::melt(tmp.data.c[,c("month",paste(tmp.outcome,"_",tmp.ration,sep = ""),"yhat")],id="month")
meltdf.c$group <- NA
meltdf.c$group[meltdf.c$variable == paste(tmp.outcome,"_",tmp.ration,sep = "")] <- "Observed values, control, "
meltdf.c$group[meltdf.c$variable =="yhat"] <- "Fitted values, control, "


meltdf.t <- reshape2::melt(tmp.data.t[,c("month",paste(tmp.outcome,"_",tmp.ration,sep = ""),"yhat")],id="month")
meltdf.t$group <- NA
meltdf.t$group[meltdf.t$variable == paste(tmp.outcome,"_",tmp.ration,sep = "")] <- "Observed values, treatment, "
meltdf.t$group[meltdf.t$variable =="yhat"] <- "Fitted values, treatment, "

finaldf <- rbind( meltdf.c,  meltdf.t)


#reshape to get se
meltdf.se.t <- reshape2::melt(tmp.data.t[,c("month",paste(tmp.outcome,"_",tmp.ration,"_se",sep = ""),"yhat_se")],id="month")
meltdf.se.t$se <- 0
meltdf.se.t$se[meltdf.se.t$variable =="yhat_se"] <- meltdf.se.t$value[meltdf.se.t$variable =="yhat_se"]  #leave the rest zero, we don't need to plot se of observed
meltdf.se.t$value <- NULL
meltdf.se.t$group <- NA
meltdf.se.t$group[meltdf.se.t$variable == paste(tmp.outcome,"_",tmp.ration,"_se",sep = "")] = "Observed values, treatment, "
meltdf.se.t$group[meltdf.se.t$variable =="yhat_se"] <- "Fitted values, treatment, "
meltdf.se.t$variable <- NULL

meltdf.se.c <- reshape2::melt(tmp.data.c[,c("month",paste(tmp.outcome,"_",tmp.ration,"_se",sep = ""),"yhat_se")],id="month")
meltdf.se.c$se = 0
meltdf.se.c$se[meltdf.se.c$variable =="yhat_se"] <- meltdf.se.c$value[meltdf.se.c$variable =="yhat_se"]  #leave the rest zero, we don't need to plot se of observed
meltdf.se.c$value <- NULL
meltdf.se.c$group <- NA
meltdf.se.c$group[meltdf.se.c$variable == paste(tmp.outcome,"_",tmp.ration,"_se",sep = "")] <- "Observed values, control, "
meltdf.se.c$group[meltdf.se.c$variable =="yhat_se"] <- "Fitted values, control, "
meltdf.se.c$variable <- NULL


finaldf.se=rbind(meltdf.se.c,  meltdf.se.t)


plot.dat1 <- merge(finaldf,finaldf.se, by = c("month","group"))
plot.dat1$lwr <- plot.dat1$value - 1.96*plot.dat1$se
plot.dat1$upr <- plot.dat1$value + 1.96*plot.dat1$se

Upper=max(plot.dat1$upr)
#######################################################################################
# Un reconciled
#######################################################################################
r <- 2
tmp.ration <- ration.value_total[r]
tmp.ration1 <- ration.name_total[r]

#treatment

tmp.data.c <- readstata13::read.dta13(paste(SurveyDataDir,tmp.ration,"_",tmp.outcome,"_meanreceipts_reconciliation_time_series_treatment0.dta",sep = ""))
tmp.data.t <- readstata13::read.dta13(paste(SurveyDataDir,tmp.ration,"_",tmp.outcome,"_meanreceipts_reconciliation_time_series_treatment1.dta",sep = ""))


#reshape the dataset to include values for month and whether predicted/observed
meltdf.c <- reshape2::melt(tmp.data.c[,c("month",paste(tmp.outcome,"_",tmp.ration,sep = ""),"yhat")],id="month")
meltdf.c$group <- NA
meltdf.c$group[meltdf.c$variable == paste(tmp.outcome,"_",tmp.ration,sep = "")] <- "Observed values, control, "
meltdf.c$group[meltdf.c$variable =="yhat"] <- "Fitted values, control, "


meltdf.t <- reshape2::melt(tmp.data.t[,c("month",paste(tmp.outcome,"_",tmp.ration,sep = ""),"yhat")],id="month")
meltdf.t$group <- NA
meltdf.t$group[meltdf.t$variable == paste(tmp.outcome,"_",tmp.ration,sep = "")] <- "Observed values, treatment, "
meltdf.t$group[meltdf.t$variable =="yhat"] <- "Fitted values, treatment, "

finaldf <- rbind( meltdf.c,  meltdf.t)


#reshape to get se
meltdf.se.t <- reshape2::melt(tmp.data.t[,c("month",paste(tmp.outcome,"_",tmp.ration,"_se",sep = ""),"yhat_se")],id="month")
meltdf.se.t$se <- 0
meltdf.se.t$se[meltdf.se.t$variable =="yhat_se"] <- meltdf.se.t$value[meltdf.se.t$variable =="yhat_se"]  #leave the rest zero, we don't need to plot se of observed
meltdf.se.t$value <- NULL
meltdf.se.t$group <- NA
meltdf.se.t$group[meltdf.se.t$variable == paste(tmp.outcome,"_",tmp.ration,"_se",sep = "")] <- "Observed values, treatment, "
meltdf.se.t$group[meltdf.se.t$variable =="yhat_se"] <- "Fitted values, treatment, "
meltdf.se.t$variable <- NULL

meltdf.se.c <- reshape2::melt(tmp.data.c[,c("month",paste(tmp.outcome,"_",tmp.ration,"_se",sep = ""),"yhat_se")],id="month")
meltdf.se.c$se <- 0
meltdf.se.c$se[meltdf.se.c$variable =="yhat_se"] <- meltdf.se.c$value[meltdf.se.c$variable =="yhat_se"]  #leave the rest zero, we don't need to plot se of observed
meltdf.se.c$value <- NULL
meltdf.se.c$group <- NA
meltdf.se.c$group[meltdf.se.c$variable == paste(tmp.outcome,"_",tmp.ration,"_se",sep = "")] <- "Observed values, control, "
meltdf.se.c$group[meltdf.se.c$variable =="yhat_se"] <- "Fitted values, control, "
meltdf.se.c$variable <- NULL


finaldf.se=rbind(meltdf.se.c,  meltdf.se.t)


plot.dat2 <- merge(finaldf,finaldf.se, by = c("month","group"))
plot.dat2$lwr <- plot.dat2$value - 1.96*plot.dat2$se
plot.dat2$upr <- plot.dat2$value + 1.96*plot.dat2$se
Lower <- min(plot.dat2$lwr)

#######################################################################################
# Combine
#######################################################################################
plot.dat2$group <- paste(plot.dat2$group,"non-reconciled")
plot.dat1$group <- paste(plot.dat1$group,"reconciled")
plot.dat <- rbind(plot.dat1,plot.dat2)

#######################################################################################
##  Plot
#######################################################################################
ggplot(plot.dat, aes(x=month, y=value, group = group, colour= group, linetype =group)) +
  geom_line(size= 1.5)+ 
  geom_ribbon(aes(ymin=lwr,ymax=upr), alpha = 0.1, colour=NA) +  
  scale_colour_manual(values=c(fColor,darkColor,fColor,darkColor,mColor,lightColor,mColor,lightColor)) +  
  scale_linetype_manual(values=c("dashed","dashed","solid","solid", "dotted","dotted","dotted","dotted"))+
  xlab("Month") +
  ylab("Average receipt value per rationcard")+ 
  guides(colour = guide_legend(override.aes = list(shape = NA))) + 
  guides(colour=guide_legend(ncol=2))+
  theme(axis.text=element_text(size=16), legend.key = element_blank(),
        axis.title=element_text(size=16),legend.position= "bottom",
        legend.title=element_blank(), legend.text=element_text(size=16), legend.spacing.x = unit(0.3, 'cm'),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank()) +
  scale_y_continuous(limit=c(Lower  ,550),breaks = seq(0,550,100)) +
  scale_x_continuous(limit=c(1,11), breaks = seq(1,11,1), labels = month.list) +
  theme(legend.key = element_blank(), legend.key.width = unit(1.5, "cm"))+
  geom_rect(data=rect_r, aes(xmin=xmin, xmax=xmax, 
                             ymin=ymin, ymax=ymax), 
            color="gray100", alpha=0.05, inherit.aes = F,
            linetype = 0)



ggsave(file=paste(OutputDir, "Figure3_PanelB.pdf", sep=""), width=12, height=8)



